home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 46 / Amiga Format CD46 (1999-10-20)(Future Publishing)(GB)[!][issue 1999-12].iso / -in_the_mag- / reader_requests / scilab / demos / velpic / acoustic.sci < prev    next >
Text File  |  1999-09-16  |  8KB  |  258 lines
  1. //[pt,dx,dz,dt]=acoustic(vel,tf,fc,spos,dx,dz,dt)
  2. //[pt[,dx,dz,dt]]=acoustic(vel,tf,fc,spos[,dx,dz][,dt])
  3. //////////////////////////////////////////////////////
  4. //                                                  //
  5. //  Scilab program to simulate forward propagation  //
  6. //  of acoustic waves with absorbing boundary       //
  7. //  conditions where p is the wavefield, s is the   //
  8. //  source, and v is the velocity field:            //
  9. //                                                  //
  10. //                      2  __ 2                     //
  11. //             P   =  v    \/   P + S               // 
  12. //              tt                                  //
  13. //                                                  //
  14. //////////////////////////////////////////////////////
  15. // vel  :Velocity distribution (matrix which represents
  16. //      :distribution of velocity as function of offset x 
  17. //      :and depth z.  Increasing offset and depth goes
  18. //      :as the increasing indices of the matrix vel).
  19. // tf   :Final time (initial time is zero)
  20. // fc   :Center frequency of wavelet (derivative of gaussian)
  21. // spos :source postion which is a 2-vector of integers 
  22. //      :where 1<=spos(1)<=nr, 1<=spos(2)<=nc 
  23. //      :and [nr,nc]=size(vel).
  24. // dx   :Sampling step in offset 
  25. // dz   :Sampling step in depth
  26. // dt   :Sampling step in time.  
  27. // ntbl :Name table containing names of all the files 
  28. //      :created containing data
  29. //
  30. //The parameters dx, dz, and dt are optional and are
  31. //calculated automatically, when not specified, to satisfy
  32. //stability conditions and to impose an acceptable level
  33. //of numerical dispersion.
  34. //
  35. //!
  36. //author: C. Bunks     date: 29-Oct-90
  37.  
  38.    [lhs,rhs]=argn(0);
  39.    lines(0);
  40.  
  41. //velocity parameters
  42.  
  43.    [nr,nc]=size(vel);
  44.    vmax=maxi(vel);vmin=mini(vel);
  45.  
  46. //default check
  47.  
  48.    if rhs=4 then,//auto calculation of dx, dz, and dt
  49.       dx=vmin/(16*fc);
  50.       dz=dx;
  51.       dmin=mini([dx,dz]);dmax=maxi([dx,dz]);   
  52.       dt=.95*dmin/(vmax*sqrt(2));//stability condition 
  53.    end,
  54.    if rhs=5 then,//auto calculation of dx and dz
  55.       dx=vmin/(16*fc);
  56.       dz=dx;
  57.       dmin=mini([dx,dz]);dmax=maxi([dx,dz]);   
  58.    end,
  59.    dmin=mini([dx,dz]);dmax=maxi([dx,dz]);   
  60.    if rhs=6 then,//auto calculation of dt
  61.       dt=.95*dmin/(vmax*sqrt(2));//stability condition 
  62.    end,
  63.  
  64. //inform user of default choices
  65.  
  66.    write(%io(2),' '),
  67.    write(%io(2),' CHOICES OF DX, DZ, AND DT:'),
  68.    print(6,[dx,dz,dt]),
  69.    write(%io(2),' '),
  70.  
  71. //ERROR CHECKING
  72. //source location
  73.  
  74.    eflag='on';
  75.    if 1<=spos(1) then, if spos(1)<=nr then,
  76.    if 1<=spos(2) then, if spos(2)<=nc then,
  77.       eflag='off';
  78.    end,end,
  79.    end,end,
  80.    if eflag='on' then,
  81.       write(%io(2),'                               '),
  82.       write(%io(2),'*************ERROR*************');
  83.       write(%io(2),'                               '),
  84.       write(%io(2),' SOURCE POSITION OUT OF BOUNDS ');
  85.       write(%io(2),'                               '),
  86.       write(%io(2),'*************ERROR*************');
  87.       write(%io(2),'                               '),
  88.       error('sim.bas'),
  89.    end,
  90.  
  91. //space discretization
  92.  
  93.    if dmax>vmin/(15*fc) then,
  94.       write(%io(2),'                                       '),
  95.       write(%io(2),'****************WARNING****************');
  96.       write(%io(2),'                                       '),
  97.       write(%io(2),' NUMERICAL DISPERSION LIKELY           '),
  98.       write(%io(2),' VARIABLES SHOULD SATISFY:             '),
  99.       write(%io(2),'                                       '),
  100.       write(%io(2),'     max([dx,dz]) < vmin/(15*fc)       '),
  101.       write(%io(2),'                                       '),
  102.       write(%io(2),'****************WARNING****************');
  103.       write(%io(2),'                                       '),
  104.    end,
  105.  
  106. //time discretization and stability check
  107. //           v*dt/dx < 1/sqrt(2)
  108.  
  109.    if (vmax*dt/dmin)>(1/sqrt(2)) then,
  110.       write(%io(2),'                                       '),
  111.       write(%io(2),'*****************ERROR*****************');
  112.       write(%io(2),'                                       '),
  113.       write(%io(2),' UNSTABLE CHOICES: vel, dt, dx, and dz '),
  114.       write(%io(2),' VARIABLES MUST SATISFY:               '),
  115.       write(%io(2),'                                       '),
  116.       write(%io(2),'    max(vel)*dt*sqrt(2) < min(dx,dz)   '),
  117.       write(%io(2),'                                       '),
  118.       write(%io(2),'*****************ERROR*****************');
  119.       write(%io(2),'                                       '),
  120.       error('sim.bas'),
  121.    end,
  122.    t=0:dt:tf;
  123.  
  124. //pre-calculation
  125.  
  126.    v2dt=(dt*dt)*(vel.*vel);
  127.    zr=0*ones(1,nc);zc=0*ones(nr,1);
  128.  
  129. //integrate forward 
  130.  
  131. //   unit=file('open','pt.dat','unknown','unformatted');
  132.    pt=[];
  133.    pkm1=0*ones(vel);pkm2=0*ones(vel);//initial conditions
  134.    for tk=t,
  135.       src=shot(tk,fc);
  136.       pk=integrate(tk,pkm1,pkm2,src,spos);
  137.       pt=[pt;pk];
  138. //      writb(unit,pk);
  139.       pkm2=pkm1;pkm1=pk;
  140.    end,
  141. //   file('rewind',unit);
  142. //   pt=readb(unit,nr*maxi(size(t)),nc);
  143. //   file('close',unit);
  144.  
  145. //end
  146.  
  147. //[utp1]=integrate(t,ut,utm1,src,spos)
  148. //[utp1]=integrate(t,ut,utm1,src,spos)
  149. //compute second order time update of acoustic wave equation
  150. //(with absorbing boundaries)
  151. // t     :current time
  152. // ut    :wavefield at time t
  153. // utm1  :wavefield at time t-dt
  154. // src   :source value at time t
  155. // spos  :source position
  156. // utp1  :wavefield at time t+dt
  157. //
  158. //!
  159. //author: C. Bunks     date: 29-OCT-90
  160.  
  161.    write(%io(2),'t='+string(t));
  162.  
  163. //calculate laplacian in the interior of the medium
  164.  
  165.    utp1= 2*ut...
  166.         -utm1...
  167.         +v2dt.*(([ut(:,2:nc) zc]+[zc ut(:,1:nc-1)]-2*ut)/dx**2...
  168.                +([ut(2:nr,:);zr]+[zr;ut(1:nr-1,:)]-2*ut)/dz**2);
  169.  
  170. //calculate boundary conditions on edges  (fix velocities)
  171. //(see paper by Reynolds, Geophysics, v. 43, 1978, pp1099-1110)
  172.  
  173. //right side boundary
  174.    ua=ut(:,nc-1);ub=ut(:,nc);uc=utm1(:,nc-2);ud=utm1(:,nc-1);
  175.    utp1(:,nc)=(ua+ub-ud)-(dt/dx)*vel(:,nc).*(ub-ua+uc-ud);
  176.  
  177. //left side boundary
  178.    ua=ut(:,2);ub=ut(:,1);uc=utm1(:,3);ud=utm1(:,2);
  179.    utp1(:,1)=(ua+ub-ud)-(dt/dx)*vel(:,1).*(ub-ua+uc-ud);
  180.  
  181. //bottom boundary
  182.    ua=ut(nr-1,:);ub=ut(nr,:);uc=utm1(nr-2,:);ud=utm1(nr-1,:);   
  183.    utp1(nr,:)=(ua+ub-ud)-(dt/dz)*vel(nr,:).*(ub-ua-ud+uc);
  184.  
  185. //top boundary (absorbing or free)
  186.    ua=ut(2,:);ub=ut(1,:);uc=utm1(3,:);ud=utm1(2,:);   
  187.    utp1(1,:)=(ua+ub-ud)-(dt/dz)*vel(1,:).*(ub-ua-ud+uc);
  188.  
  189. //calculate boundary conditions at corners
  190.    utp1(1,1)=ut(1,1)+((dt/dx)*vel(1,1))*(ut(2,2)-ut(1,1));
  191.    utp1(1,nc)=ut(1,nc)+((dt/dx)*vel(1,nc))*(ut(2,nc-1)-ut(1,nc));
  192.    utp1(nr,1)=ut(nr,1)+((dt/dx)*vel(nr,1))*(ut(nr-1,2)-ut(nr,1));
  193.    utp1(nr,nc)=ut(nr,nc)+((dt/dx)*vel(nr,nc))*(ut(nr-1,nc-1)-ut(nr,nc));
  194.  
  195. //add in source
  196.    sz=spos(1);sx=spos(2);
  197.    utp1(sz,sx)=utp1(sz,sx)+dt**2*src;
  198.  
  199. //end
  200.  
  201. //[dg]=shot(t,fc)
  202. //[dg]=shot(t,fc)
  203. //calculate shot values as a function of time
  204. //as the derivative of a gaussian:
  205. //
  206. //   g=exp(-(t-m)**2/(2*sigma**2))
  207. //   dg=-g*(t-m)/sigma**2
  208. //
  209. //Center frequency for the derivative of the gaussian is
  210. //at fc=1/(2*m). The ratio of the minimum velocity (vm) 
  211. //to center frequency (fc) (i.e., the smallest spatial wavelength)
  212. //should be between 10 and 20 times greater than the largest 
  213. //spatial discretization.  Here we take the factor to be 15:
  214. //
  215. //              vm/fc = 2*m*vm > 15*maxi(dx,dz)
  216. //or
  217. //                  m = 7.5*maxi(dx,dz)/mini(vel)
  218. //
  219. // t  :current time
  220. // fc :center frequency of the wavelet
  221. // dg :derivative of a gaussian at time t
  222. //
  223. //!
  224. //author: C. Bunks     date: 29-OCT-90
  225.  
  226.    m=1/(2*fc);
  227.    sig=m/4;
  228.    sig2=sig**2;
  229.    g=exp(-(t-m)**2/(2*sig2))/sqrt(2*%pi*sig2);
  230.    dg=-g*(t-m)/sig2;
  231.  
  232. //end
  233.  
  234. //[pt]=get_data(ntbl,entry)
  235. //Search for a data file written on disk
  236. // ntbl  :table of file names (first two entries give 
  237. //       :data dimensions)
  238. // entry :integer giving the entry in the table-2
  239. // pt    :returned data file
  240.  
  241.    ts=maxi(size(ntbl));
  242.    nr=evstr(ntbl(1));
  243.    nc=evstr(ntbl(2));
  244.    if entry<=ts-2 then,
  245.       filename=ntbl(entry+2);
  246.       unit=file('open',filename,'unknown');
  247.       pt=read(unit,nr,nc);
  248.       file('close',unit);
  249.       plot2d(pt,'agc',[0,%pi/4],'x');
  250.    else,
  251.       write(%io(2),' '),
  252.       write(%io(2),' Table Entry to Large...Max Value='+string(ts-2)),
  253.       write(%io(2),' '),
  254.    end,
  255.  
  256. //end
  257.  
  258.